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Abstract 

We investigate the extension of the multilevel Monte Carlo path 
simulation method to jump-diffusion SDEs. We consider models with 
finite rate activity , using a jump-adapted discretisation in which the 
jump times are computed and added to the standard uniform dis- 
cretisation times. The key component in multilevel analysis is the 
calculation of an expected payoff difference between a coarse path 
simulation and a fine path simulation with twice as many timesteps. 
If the Poisson jump rate is constant, the jump times are the same on 
both paths and the multilevel extension is relatively straightforward, 
but the implementation is more complex in the case of state-dependent 
jump rates for which the jump times naturally differ. 



1 Introduction 

In the Black-Scholes Model, the price of an option is given by the expected 
value of a payoff depending upon the solution of a stochastic differential 
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equation(SDE) satisfied by the stock price. The model assumes that the 
behavior of the stock price is depicted by a SDE driven by Brownian motion, 

dS(t)=a(S,t)dt + b(S,t)dW(t), < t < T, (1) 

with given initial data S . 

Although this model is widely used, the fact that asset returns are not 
log-normal has motivated people to suggest models which better capture the 
characteristics of the stock price dynamics. Merton pVIer76] proposed a jump- 
diffusion process for the stock price. To be specific, the stock price follows a 
jump-diffusion SDE: 

dS(t) = a(S(t-),t)dt + b(S(t-),t)dW(t) + c(S(t-),t)dJ(t), < t < T, 

(2) 

where the jump term J{t) is a compound Poisson process Y^i=i \Xi ~ -Q> the 
jump magnitude Yj has a prescribed distribution, and N(t) is a Poisson 
process with intensity A, independent of the Brownian motion. Due to the 
existence of jumps, the process is a cadlag process, i.e. having right continuity 
with left limits. We note that S(t—) denotes the left limit of the process 
while S(t) = lim s ^ t+ S(t). In [Mer76j, Merton also assumed that logYJ has 
a normal distribution with mean a and variance b, namely log Yj ~ N(a, b). 

There are several ways to generalize Merton model from different aspects. 
A possible way is to consider the case where the frequency of jump is infinite, 
where general Levy processes can be used. Another direction (for example 
in |GM03| ) is to introduce dependency between parameters, leaving the ex- 
pected number of jumps finite within finite horizon. As a particular numerical 
examples, we consider the case where instantaneous jump rate relies on the 
stock price, namely \{t) = X(S t ,t). 

In pursuit of risk-neutral pricing of options, we are interested in the ex- 
pected value of a function of the terminal state, f(S(T)). In the simple case 
of a European option, the expectation can be directly simulated, while in 
the case of Asian, lookback and barrier options the valuation depends on 
the entire path S(t),0 <t <T. The expected value can be estimated by a 
simple Monte Carlo method with a proper numerical discretisation scheme. 
However, to achieve a root-mean-square (RMS) error of 0(e) using an Euler- 
Maruyama discretisation would require 0(e~ 2 ) independent paths, each with 
0(e _1 ) timesteps, leading to a computational complexity of 0(e -3 ). This is 
quite time consuming compared to the case of path-independent options. 

Giles [Gil07j |Gil08bj has recently introduced a multilevel Monte Carlo 
path simulation method for the option pricing calculation. This improves 
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the computational efficiency of Monte Carlo path simulation by combining 
results using different numbers of timesteps. This can be viewed as a gen- 
eralisation of the two-level method of Kebaier |Keb05] and is also similar in 
approach to Heinrich's multilevel method for parametric integration |Hei01j . 
The first paper [Gil08bJ proposed the multilevel Monte Carlo approach and 
proved that it can lower the computational complexity of path-dependent 
Monte Carlo evaluations to 0(e _2 (loge) 2 ), verified by numerical results us- 
ing the simple Euler-Maruyama discretisation. The second paper |Gil07j 
demonstrated that the computational cost can be further reduced to 0(e~ 2 ) 
by using the Milstein discretisation. This has been extended by Dereich and 
Heidenreich |DH11[ IDerllj to approximation methods for both finite and in- 
finite activity Levy-driven SDEs with globally Lipschitz payoffs. The work in 
this paper differs in considering simpler finite activity jump-diffusion models, 
but more challenging non-Lipschitz payoffs, and also uses a more accurate 
Milstein discretisation to achieve an improved order of convergence for the 
multilevel correction variance which will be defined later. 

In this paper we applies the multilevel approach to the Monte Carlo simu- 
lation of path-dependent option pricing under jump-diffusion processes. We 
first consider the case where the jump rate is constant then take into ac- 
count the state-dependent rate case. In both cases, in order to calculate 
coarse-path samples from fine-path sample using brownian interpolation , we 
adopt a jump-adapted Milstein discretisation scheme proposed by [Pla82]. 
which explicitly simulates the times when jumps occur. Furthermore, we con- 
struct multilevel estimators for corresponding path-dependent payoffs coping 
with challenges caused by jumps. Through constructing payoff estimators 
by Brownian bridge technique, high order multilevel correction term vari- 
ance convergence rate is achieved. In the state-dependent rate case, we use 
two approaches, which are called cumulative intensity method and thinning 
method to tackle the unsynchronization of jump times in the fine and coarse 
grids. Numerical results show similar improvement in computational effi- 
ciency compared with previous achievement |Gil07j for diffusion processes. 
Generally, using the jump-adapted Milstein scheme with the multilevel ap- 
proach, we can reduce the computation cost to 0(e~ 2 ) in terms of RMS error 
e. 

In the following parts of the paper, we first review the Multilevel Monte 
Carlo method for diffusion processes. The next section describes the jump- 
adapted discretisation of jump-diffusion processes and its advantages for fa- 
cilitating the multilevel approach. Then we discuss the path simulation and 
estimator construction for the jump-adapted discretisation with the mul- 
tilevel approach and present numerical results of Asian, lookback, barrier 



3 



and digital options. The next part establishes two methods to deal with 
state-dependent intensity. The final section draws conclusions and indicates 
directions of future research. 



2 Multilevel Monte Carlo method 

Suppose we perform Monte Carlo path simulations with fixed grid timesteps 
he = 2~ e T, I = 0, 1, . . . , L. For a given Brownian path W(t), let P denote 
the payoff, and let Pe denote its approximation by a numerical scheme with 
timestep he- As a result of the linearity of the expectation operator, 

L 

E[P L ] = E[P ] +X>[^-^-il- (3) 

t=\ 

Let Y denote an estimator for E[P ] using N paths. Suppose for different 
/ > 0, we use Nt independent paths to estimate E[P^ — p£_ x ]. 

% = N^(P®-P® 1 ). (4) 

i=i 

The Multilevel method facilitates the fact that V[P^ — Pg-i] decreases with I 
to adaptively choose N( and hence reduce the computational cost. 

The cost reduction effect is summarized in the following theorem: 

Theorem 2.0.1. Let P denote a functional of the solution of stochastic dif- 
ferential equation (T7]j for a given Brownian path W(t), and let Pe denote the 
corresponding approximation using a numerical discretisation with timestep 
h e = 2~ l T. 

If there exist independent estimators Ye based on Ne Monte Carlo samples, 
and positive constants a> |, j3, ci, C2, C3 such that 



1) E[P e -P) 
u) We] = 



< ci hi 

E[P ], 1 = 

E[Pe-Pe-il l>0 

Hi) W[Y e ] < c 2 N- x h{ 
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iv) Ci, the computational complexity ofYg, is bounded by 

C t < c 3 N e hj\ 



then there exists a positive constant c 4 such that for any e<e 1 there are 
values L and Nt for which the multilevel estimator 



£=0 



has a mean-square-error with bound 



MSE = E 



Y - E[P] 



< 6 



with a computational complexity C with bound 

c 4 e" 2 , P > 1, 

C< { c 4 e- 2 (log e ) 2 , = 1, 

C4e -2- M /a ) < /3 < 1. 



Proof. See |Gil08bj . 



□ 



In the case of the jump-adapted discretisation, he should be taken to 
be the uniform timestep at level I, to which the jump times are added to 
form the set of discretisation times. We have to define the computational 
complexity as the expected computational cost since different paths may have 
different numbers of jumps. However, the expected number of jumps is finite 
and therefore the cost bound in assumption iv) will still remain valid for an 
appropriate choice of the constant c 3 . 

According to the theorem, the larger the variance convergence rate (3 , the 
greater the reduction is the computation cost by the multilevel algorithm. 
In the case of a Lipschitz continuous European payoff, using the Milstein 
discretisation immediately leads to the result that Ve = O(hj), corresponding 
to (3 = 2. Thus the main task to improve the performance of the multilevel 
method is to use using schemes with high order strong convergence rate and 
constructing appropriate estimators so that > 1 could be achieved. For the 
jump-diffusion process, the objective is to obtain (3 > 1 through adopting a 
high order scheme. 
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3 A Jump-adapted Milstein discretisation 



To simulate jump-diffusion processes, it is possible to use fixed time grid 
schemes as for geometric Brownian motion. The Euler-Maruyama scheme for 
jump-diffusion processes has 0(Vh) strong convergence ( |PlalOj ). However, 
it would be more difficult to pursue higher order strong convergence. To 
achieve a higher order strong convergence for jump- diffusion processes, the 
It 6- Taylor expansion will involve some double integrals of white noise and 
the Poisson random measure |BLP05] , which increases the complexity of the 
simulation. 

Another problem which might be encountered for fixed-time grid schemes 
is the construction of estimators for the payoff function of path-dependent 
options. Adoption of the previous Brownian bridge interpolation is difficult 
since the minimum or other functional of paths is difficult to calculate since 
the joint density of diffusion and jump is much more complex than pure 
diffusion one. 

In order to avoid simulating double stochastic integrals as well as to iden- 
tify the time at which the jump occurs, we use the so-called jump-adapted 
approximation proposed by Platen in |Pla82] . This jump-adapted scheme 
would improve the computational tractability compared to other fixed time 
grid discretisation schemes with the same weak/strong convergence order. 

Suppose that we have simulated the jump time grid J = {ti, r 2 , . . . , r m }, 
which includes times at which jumps occur in the [0, T\. On the other hand, 
consider a fixed time grid constituted by iV timesteps, = i x |, % — 
1, . .. , N, which is used in discretisation schemes of Brownian SDEs. Now 
consider a superposition of them as a new grid T = {0 = t < t\ < t 2 < 
. . . < t M = T}. As a result, the length of timestep of the new grid will be 
no greater than h = -S. 

Within every timestep of the new grid, the diffusion part is separated 
from the effect of the possible jump, because the jump only occurs at the 
grid point. Thus we can approximate the path with established schemes 
for diffusion processes, and deal with corresponding adjustment if there is a 
jump at the right end of interval. This procedure is called the jump-adapted 
scheme. 

The algorithm of simulation via the jump-adapted scheme could be de- 
scribed as the following steps: 

1. Set i = l, j = 1, if = t' ; 
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2. Generate jump time Tj in terms of its distribution; 

3. While (n < t'j) do 

(1) . Simulate the process within \t',Tj), in which the process is driven 
purely by Brownian motion, then simulate the jump at Tj; 

(2) . Set the timestep hi = t { — t', t' = m 

(3) . % — and generate next jump time Tj in terms of its distribution; 

4. Simulate the process within [t',^]; 

5. Set the timestep hi — Tj — t' , j — j + 1, t' — and goto 3. 

Now we introduce a jump-adapted Milstein scheme for a scalar jump- 
diffusion SDE 



S n + a n h n + b n AW n + \-^b n (AW* - h n ), 

| 5+ 1 + c ^"+l> tn+i){Yi - 1), when t n+x = r,; 
1 S~ +1 , otherwise. 

Where the subscript n is used to denotes the timestep index, S~ = S(t n —) 
is the left limit of the approximated solution, and Yi is the jump magnitude 
at Tj . 

In sum, jump-adapted schemes explicitly compute jump times, which are 
relatively rare in the entire time span. Thus, compared to fixed time schemes, 
they save the computation cost for generating Poisson random numbers when 
the timestep tends to zero. Furthermore, as we will see later, in terms of path 
simulation, jump- adapted discretisations have a very crucial property that 
keeps the multilevel approach valid: within each timestep we can neglect the 
jump term and only take the Brownian component into consideration. As 
a matter of fact, the scheme can conveniently adopt the Brownian bridge 
technique used for estimator construction in the previous paper |Gil07j so 
that improved convergence could be obtained as well. 

4 Multilevel approach in the presence of jump 

In all of the cases to be presented, we simulate the paths using the jump- 
adapted Milstein scheme proposed above. 



S 



n+l 



S n 



+1 
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In the case of a jump-diffusion process, since theorem of computation 
complexity in the case of diffusion processes requires the weak convergence 
and ML estimator convergence of discretisation schemes, we have to justify 
them accordingly for different construction of estimators. These numerical 
analysis is being done in a working paper in preparation |XG11] . 

Apart from theoretical issues, there arise two challenges in the imple- 
mentation. The first problem is that the path simulation on coarse levels 
needs to be revised in the presence of varying timesteps of the jump-adapted 
time grid. Another is how to devise suitable estimators for various payoffs in 
coping with a jump- adapted time grid. The third concern is whether opti- 
mal samples on each level used by previous algorithm in |Gil08b] should be 
modified. Due to presence of jump, the computational cost is 

L N e 

cost = j2J2( N T+ 2t )> 

e=o i=i 

where is the number of jumps in each scenario. However, the expected 
number of jumps is finite and therefore the cost bound in assumption iv) will 
still remain valid for an appropriate choice of the constant C3 therefore using 
previous algorithm should work well. In implementation numerical results 
indicates that this is appropierate for small A. 

4.1 Path simulation of multilevel approach in the jump- 
adapted time grid 

When the multilevel approach is applied to the simulation with a fixed time 
discretisation of a jump-diffusion process, the algorithm maintains the orini- 
gal framework straightforward. While for jump-adapted schemes, construc- 
tion of coarse path simulation needed for the estimators (j3j) of the payoff 
differs from previous case. 

In the case of fixed time grid for geometric Brownian motion, every path 
sample used for calculating Pg — P^_i comes from two discrete approxima- 
tions of Brownian path, which are called the fine path and the coarse path. 
For every I, every timestep of coarse grid is completely the same as the cor- 
responding two timesteps of the fine grid starting from the same endpoint. 
Since the brownian increment of the coarse timestep is equidistributed to the 
sum of two increments of corresponding fine ones, we can do the coarse path 
simulation without generating extra random normal number. 

While in the case of jump-adapted time grid, due to the presence of 
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Figure 1: Midpoint construction 



jump time, the construction of a path sample in the coarse grid using the 
brownian increments generated in the fine grid needs to be clarified. In 
this case, the path sample is discontinuous in its jump time. To avoid such 
discontinuity, notice that within each timestep of jump-adapted time grid, 
the path is purely driven by Brownian component and therefore reserves 
continuity. For the coarse grid and the fine grid in the same level, call the 
finer grid points midpoints for short. For a particular midpoint, the timestep, 
which is formed by the last jump time before this midpoint and the first 
jump time after it is called midpoint timestep. In this timestep, Brownian 
increment of the coarse path sample can be obtained by the summation of 
Brownian increments of corresponding two timestep in fine grid. In remaining 
timesteps of coarse grid, construction of the coarse path sample uses the same 
Brownian increment as the one in fine grid does. Hence we have defined the 
coarse path sample construction according to this midpoint construction, 
which can be seen clearly in figure [TJ 



4.2 Estimator construction 

For the estimator construction, there comes concern whether extra bias is 
introduced for estimator in the jump-diffusion process. To secure the cor- 
rectness of the identity we must avoid introducing any undesired bias, hence 
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it is required that 



E[P/] = E[P?\. 



(5) 



This means that the definitions of Pi when estimating E\Pi — Pe-i] and 
E[p£+i — Pi] must have the same expectation. 

In the case of path-dependent payoffs in geometric Brownian motion, we 
approximate the payoff function by Brownian bridge interpolation technique 
using path values on the fine and the coarse grid. This technique is also 
available in jump case to reduce the variance of the estimator, where the 
coarse-path estimator will involve the information from generations of fine- 
path sample in the two corresponding timesteps. One thing to notice in jump- 
adapted schemes is that we can utilise this construction only for the timesteps 
including midpoint. For other timesteps of coarse grid, the construction of 
coarse-grid estimators will be the same as the fine-grid one. Estimators will 
be discussed in the following respectively corresponding to their payoffs. 

In the following we will show the numerical results of several options. All 
of them are done for Merton model in which the jump-diffusion SDE under 
risk-neutral measure is 



where A is the jump intensity and jump magnitude satisfies log Yj ~ N(a, b), 
r is the risk-free rate, a is the volatility of stock price and m = E[Yj] — 1 is 
the compensator to ensure the discounted stock price is a martingale. All of 
the simulations in this section use the parameter values So = 100, K — 100, 
T = 1, r = 0.05, a = 0.2, a = 0.1, b = 0.2, A = 1. We thank Giles providing 
the code for [Gil08aj . based on which we can produce the current code to 
generate numerical results and figures. 

4.3 Vanilla call option 

For the vanilla option with the payoff exp(— rT) max.(S(T) — K, 0), Figure [2] 
shows the numerical results. 

The top left plot shows the behaviour of the variance of both Pi and the 
multilevel correction Pi—Pi^i, estimated using 10 5 samples so that the Monte 
Carlo sampling error is negligible. The slope of the MLMC line indicates that 
Vi = Y[Pi— Pe-i] =0(hj), corresponding to (3 = 2 in condition in) of Theorem 
12.0.11 The top right plot shows that K[Pi — Pi-i] is approximately O(hi), 
corresponding to a = 1 in condition i). Noting that the payoff is Lipschitz, 



dS(t) 
W 7 ) 



( y r-Xm)dt + adW(t) + dJ(t), < t < T. 
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Figure 2: Vanilla option 



both of these are consistent with the first order strong convergence proved 
in jPlalOj . 

The bottom two plots correspond to five different multilevel calculations 
with different user-specified accuracies to be achieved. These use the numer- 
ical algorithm given in [Gil08bJ to determine the number of grid levels, and 
the optimal number of samples on each level, which are required to achieve 
the desired accuracy. We use the computational cost Yl^i=oYl!ili{^T > + ^) 
to take account into the effect of jump. The left plot shows that in each case 
many more samples are used on level than on any other level, with very few 
samples used on the finest level of resolution. The right plot shows that the 
the multilevel cost is approximately proportional to e -2 , which agrees with 
the computational complexity bound in Theorem 12.0.11 for the /3 > 1 case. 
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4.4 Asian option 



The payoff of the Asian option we consider is 

P = exp(— rT) max (0, S—K) , 

where 

~S = T _1 I S(t) dt. 
Jo 

rirp = T/h is the number of timesteps. |Gil07] shows that accuracy can be 
achieved by approximating the behaviour of a process within a timestep as 
an Ito process with constant drift and volatility, conditional on the computed 
endpoint values S n . Taking b n to be the constant volatility within the interval 
[t n ,t n+ i], in other words, we define brownian interpolation in the coarse grid 
at t as 

S(t) = S n + fi(S~ +1 - S n ) + b n [W t -W n - fi(W n +i - W n )l (6) 
where // = (£- t n )/h, h = t n+1 - t n . 
This implies 

rtn+i 

/ S(t) dt = \h{S{t n ) + S(t n+1 -)) + b n AI n , 

where AI n is 

l-tn + l 

AI n := / (W(t) - W{t n )) dt - | hAW, 

satisfying AI n ~ N(0, h 3 /12) , and is independent of AW. Let b n = b(S n , t n ), 
the fine-path approximated payoff would be 

ray — 1 

S f = T- 1 ( \ h (Sn+S; i+l ) + b n All) . 

ra=0 

In a jump-adapted grid, the coarse-path approximation is the same in 
most timesteps except in the midpoint timestep AI° is derived from the 
fine-path values, namely 
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Figure 3: Asian option 



rtn+2 

/ (W(t) - W(t n )) dt - \{t n+2 - t n )(W(t n+2 ) - W(Q) 

Jt n 

= [ " +1 (w(t) - W(t n )) dt - \ (t n+1 - t n ) (W(t n+l ) - W{t n )) 

pt n +2 

+ / (W(t) - W(t n+1 )) dt - \ (t n+2 - t n+1 ) (W(t n+2 ) - W(t n+1 )) 

Jt n +h 

+ \ (tn+2 - tn+l) (W(t n+1 ) - W(t n )) - \ (t n+1 - t n ) (W(t n+2 ) - W(t n+1 )) , 

and thus 

AF = All 1 + A/f + \ (t n+2 - t^AW? 1 - (1 - ^AW^ 2 ), 

where /z n = (t n+2 — t n+1 )/(t n+2 — t n ), AI C is the value for the coarse timestep 
in midpoint timestep; AP 1 and AW^ 1 are the values for fine timestep in the 
first fine-path timestep constituting the midpoint timestep; AP 2 and AW f2 
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are the values for the fine timestep in the latter one constituting the midpoint 
timestep. 

Figure [3] shows the numerical results for parameters 5(0) = 100, i^ = 100, 
T— 1, r = 0.05, a = 0.2, a = 0.1, b = 0.2, A = 1. All the results are similar to 
the pure diffusion case. 



4.5 Lookback option 

The payoff of the lookback option we consider is 



P = exp(-rT) ^(T)-mm5(i) 

Previous work [Gil07] achieved a second order convergence rate for the mul- 
tilevel correction variance using the Milstein discretisation and an estimator 
constructed by approximating the behaviour within a timestep as an Ito pro- 
cess with constant drift and volatility, conditional on the endpoint values 
S n and S n+ i. Brownian Bridge results (see section 6.4 in [Gla04]) give the 
minimum value within the timestep [t n ,t n+ i], conditional on the end values, 
as 



S n ,min = \ \S n + S n+1 - y (j3 n+1 -S n ^) - 2 b 2 n h log U n \ , (7) 

where b n is the constant volatility and U n is a uniform random variable on 
[0, 1]. The same treatment can be used for the jump-adapted discretisation 
in this paper, except that S~ +1 must be used in place of S n+ i in (J7J). 

Equation (J7J) is used for the fine path approximation, but a different 
treatment is used for the coarse path, as in |Gil07] . This involves a change 
to the original telescoping sum in fl3]) which now becomes 



E[Pi]=E[pf]+j2m f -pt 



where P[ is the approximation on level I when it is the finer of the two levels 
being considered, and P$ is the approximation when it is the coarser of the 
two. This modified telescoping sum remains valid provided E[P/] = E[PJ]. 

Considering a particular timestep in the coarse path construction, we 
have two possible situations. If it does not contain one of the fine path 
discretisation times, and therefore corresponds exactly to one of the fine 
path timesteps, then it is treated in the same way as the fine path, using 
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Figure 4: Lookback option 



the same uniform random number U n . This leads naturally to a very small 
difference in the respective minima for the two paths. 

The more complicated case is the one in which the coarse timestep con- 
tains one of the fine path discretisation times t' , and so corresponds to the 
union of two fine path timesteps. In this case, the value at time t' is given 
by the conditional Brownian interpolant 

S(t') =S n + n (S- +l - S n ) + b n (W(f) -W n -n (W n+1 - W n )) , (9) 

where \i = {t' — t n )/(t n+ i— t n ) and the value of W(t') comes from the fine path 
simulation. Given this value for S(t'), the minimum values for S(t) within 
the two intervals [t n ,t'] and [t',t„ + i] can be simulated in the same way as 
before, using the same uniform random numbers as the two fine timesteps. 

The equality E[P/] = E[P £ C ] is respected in this treatment because W(t') 
comes from the correct distribution, conditional on W n+ i, W n , and therefore, 
conditional on the values of the Brownian path at the set of coarse discretisa- 
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tion points, the computed value for the coarse path minimum has exactly the 
same distribution as it would have if the fine path algorithm were applied. 

Further discussion and analysis of this is given in |XGllj . including a 
proof that the strong error between the analytic path and the conditional 
interpolation approximation is at worst 0(h log/i). 

Figure H] presents the numerical results. The results are very similar to 
those obtained by Giles for geometric Brownian motion |Gil07] . The top 
two plots indicate second order variance convergence rate and first order 
weak convergence, both of which are consistent with the 0(h log/i) strong 
convergence. The computational cost of the multilevel method is therefore 
proportional to e~ 2 , as shown in the bottom right plot. 



4.6 Barrier option 



We consider a down-and-out call barrier option for which the discounted 
payoff is 

P = exp(-rT) (S(T)-K)+ 1 Wt>b} , 

where M T = mm < t < T S(t). The jump-adapted Milstein discretisation with 
the Brownian interpolation gives the approximation 



exp(-rT) (S{T)-K) + 1 



{m t >b} 



where = min <t<T S(t). This could be simulated in exactly the same 
way as the lookback option, but in this case the payoff is a discontinuous 
function of the minimum Mt and an 0(h) error in approximating Mt would 
lead to an 0(h) variance for the multilevel correction. 

Instead, following the approach of Cont & Tankov (see page 177 in [CT04j ). 
it is better to use the expected value conditional on the values of the discrete 
Brownian increments and the jump times and magnitudes, all of which may 
be represented collectively as T . This yields 



E 



exp(-rT) (S(T)-KY\ {Tlr>B} 



E 
E 



exp(-rT) (S(T)-K) + E [l {X r T>B} 

tit— 1 

exp(-rT) (S(T) - K)+ J] % 

n=0 



T 



where p n denotes the conditional probability that the path does not cross the 
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Figure 5: Barrier option 

barrier B during the n th timestep: 

- , | -2(g„-B) + (§„- +1 -B) + . 

p " = 1 " exp| vk<u»-u) (10) 



n ) 



For fine-path value, we compute p{ where b n is defined equal to b(S n ,t 
within each timestep of the jump-adapted grid. Note that S n +i in ffTOj) should 
be replaced by the value of left limit of the endpoint S~ +l , as in the lookback 
calculation. 

For the coarse path calculation, we again deal separately with two cases. 
When the coarse timestep does not include a fine path time, then we again 
use (imp . In the other case, when it includes a fine path time t' we evaluate 
the Brownian interpolant at t' and then use the conditional expectation to 
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obtain 



Pn = < 1 - exp 



-2 (S n -B) + (S(t)-B) 

K (t'-t n ) 

-2(S(t'+(S- +1 -B 




Figure shows the numerical results for K = 100, B = 85. The top 
left plot shows that the multilevel variance is O(h^) for (5 rj 3/2. This is 
similar to the behavior for a diffusion process |Gil07j . The bottom right plot 
shows that the computational cost of the multilevel method is again almost 
perfectly proportional to e~ 2 . 



4.7 Digital option 

The digital option considered here has the discounted payoff 

P = exp(-rT) 1 {S{T)>K} . 

In |Gil07] . a multilevel variance convergence rate of O(h^ ) is achieved 
by smoothing the payoff using conditional expectation given the brownian 
increments terminated one timestep before reaching the terminal time T. 
The estimator is the probability that S nt > K under assumption of simple 
Brownian motion with constant drift a n£ _i and volatility b ni _i within last 
timestep where denotes number of fine-path timesteps: 



E[P e - P|_J = E[ E[ /(££_!) - | AW h i = 1, . . . , m - 1] ] 

where $ is the cumulative density function of Normal variable, h is fine-path 
fixed-time timestep. 

In the jump-adapted time grid, the relation between last jump time and 
the last timestep before expiry leads to different expressions of above condi- 
tional expectation estimator. Let In fact, there would be three cases: 

1. The last jump time J happens before penultimate fixed-time timestep, 
i.e. J < (N — 2)^-, where iV is the number of timesteps in fixed-time 
fine grid; 
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case 1 




case 2 




case 3 



t 

J 



Figure 6: Construction of conditional expectation estimator for Digital op- 
tion 

2. The last jump time is within the last fixed-time timestep , i.e. J > 

3. The last jump time is within penultimate fixed-time timestep, i.e. (N — 
1)5 > J>(N-2)%. 

Correspondingly, different fine-path and coarse-path estimators are shown 
in the following. 

1. In case 1, the fine-path grid and coarse-path grid is the same as the 
previous diffusion case, hence we do not need to change the estimator, 
which is the probability that S nT > K under assumption of approxi- 
mation the dynamics as a constant drift a{ T _ 1 =a(S'£ r _ 1 , T— h nT ) and 
volatility b( lT _ 1 = b(Sl T _ 11 T—hi) Brownian motion within last timestep. 
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where $ is the cumulative Normal distribution, h = T/2 e is fine-path 
fixed-time timestep. 



S c nT _ 2 +2a c nT _ 2 h+b c nT _ 2 AW h - K 



tit—2 



2. In case 2, last timestep of fine path would be hj = T — t riT _i. Due to 
discontinuity of path before last jump, we must use the same estimator 
for both fine and coarse path 



pf 



*Sn T -l + a n T -l^j K 



u n T -l 



pc (T. I n T — l n T — l J 



tit — 1 1 



3. In the last case, J denotes the last jump time, and hj = T — J — h. We 
again utilise Brownian increment generated for fine path Wh ■ 



P = $ 



b f 

u n T -l 



\b C nT - 2 \ 



K 



In three cases, the conditional expectation of coarse-path estimator is 
equal to fine-path one, thus equality ()5]) is justified. Figure |6] clearly demon- 
strates three cases. 

Figure [7] shows the numerical results for parameters S'(0) = 100, i^ = 100, 
T = 1, r = 0.05, cr = 0.2. The top left plot shows that the variance is 
approximately 0(h 3 e ), corresponding to f3 = 1.5. The reason for this is 
similar to the argument in |Gil07j . 
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Figure 7: Digital option 



A different feature compared to the geometric Brownian motion case is 
that the variance of the level estimator is a constant increasing with jump 
rate A, instead of zero. The reason is simply because that the trajectories in 
level do vary in each simulation. 



5 Path-dependent rate cases 



In the case of path-dependent jump rates, which means the jump intensity 
depends on the process, for instance A = X(S t ,t), the implementation of 
multilevel becomes difficult due to the fact that path may jump at different 
time in the fine grid and coarse grid. These differences in path sample might 
enlarge the difference of quantities used in the computation of payoff between 
the fine grid and coarse grid, such as final payoff in vanilla and digital option, 
minimum estimator by Brownian bridge interpolation in lookback option and 
crossing probability in barrier option. This leads to a increased variance of 
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every single splitted multilevel estimator and finally decrease the reduction 
of computation cost. 

To tackle this obstacle, we propose two approaches. The first one facil- 
itates the idea of change of measure in dealing with multilevel method for 
discontinuous payoffs, as we used before in dealing with digital option. The 
second one uses thinning technique to simulate the desired jump time with 
inhomogeneous rate through acceptance-rejection procedure. We again need 
to change the measure when evaluating Pi — P^_i in both grids to reduce 
variance. 



5.1 Cumulative intensity method 

In the first approach, which we call cumulative intensity method, we generate 
jump times from cumulative intensity which is computed in accordance with 
the dynamics evolution. We can see the idea from the case of deterministic 
time-inhomogeneous rate. For a time-inhomogeneous Poisson process whose 
instaneous intensity is A(s), the distribution of next jump time Tj+i given Tj 
is 

P(r i+ i -Ti< t\n) = 1 - exp(- jf ' A(s)ds). (12) 

construct jump-adapted schemes by approximating the instaneous jump 
rate adaptively with the evolution of the path sample at each time grid point. 
In other words, since jump intensity is path-dependent, the jump-adapted 
time discretisation grid should be generated corresponding to the evolution 
of the underlying process. In 

In other words, 

n+i 

A(s)ds ~ (Exponential distribution with parameter 1.) 

i 

We can use this property to generate r i+1 : 

Ar i+l = inf{t > n : J X(s)ds > S i+1 }, Si ~ 8(1). 

So 

n+i =n + Ar j+ i. 

If the cumulative intensity 

A(t) = [ \( s )ds 
Jo 
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does not admit an explicit expression or it is computational intensive, we can 
use simple quadrature with linear interpolation to generate the approximated 

n+i : %. 



X(s)ds PS ^ ^i h i := A " 



where t n is the next nth point in the fixed-timestep grid, hi is the timestep of 
the grid, Aj = A(rj +X)}=i hj) * s the instaneous intensity on the left endpoint 
of each timestep. 

Let N = inf{n : A n > then 

A %i = \ (13) 

An 

is a approximation of Arj+i. 

In the case of random intensity where A only depends on current state, 

we can define Tn ELS 



Ti := inf {t : / A s ds > E\ + • — \- £{}. 
Jo 



We can use approximation of integral and ffTBj) to generate %. 

The algorithm of jump- adapted (Milstein) scheme with cumulative inten- 
sity would be: 

Algorithm (jump-adapted Milstein scheme with cumulative intensity) 

Suppose that we have a fixed time grid constituted of N timesteps, t[ = 
i x J, i= l,...,N. 

1. Let A,t,E = 0, i, j = 1. Draw an exponential r.v. £j with parameter 

i; 

2. While t < T, do 

(a) A! = A + X(3(t),t)(f i -t),E = S j + E, 
i. If A' > Sj, 

Let ft = 1~ A = E. 
^ \{S{t),t) 

Mark Tj = t + h as a jump time. Let j = j + 1 and generate 
P.- 
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ii. Otherwise A = A' . h = t[ — t; i = i + 1. 

(b) (Use Milstein scheme) Simulate the evolution of the process from 
t to t + h, obtaining value of S(t + h) depending on whether t + h 
is a jump time. 
Let t = t + h. 

As a side output, after the finishing of algorithm we have got the ap- 
proximated jump times Tj, which combines ^ = t x |, i = 1, . . . ,N forming 
jump-adapted grid T = {0 = t < t 1 < t 2 < ■ ■ ■ < tu = T}. 

The point process N t corresponding to the stopping times is defined by 

oo 

N t := J2 %<*}• ( 14 ) 
i=i 

This process is indeed a point process with piecewise constant intensity 
\(S(t k ),t k ), k = 0,...,M-l in [t k ,t k+1 ). 

5.1.1 Multilevel treatment 

When it comes to multilevel approach, the problem is that the current in- 
tensity may be different in fine and coarse grid, which leads to different 
distributions of next jump time. This causes two problems. First, we can no 
more use the random number generated for the path increment of fine grid 
to the one of the coarse grid, which is intrinsically unacceptable for multi- 
level approach. Secondly, the different final jumps may make a big difference 
between the payoffs in the fine grid and in the coarse grid, which is a similar 
challenge for digital option. 

To handle these problems, we change the measure of Poisson rate in 
calculating the expectation in the coarse grid so that the distribution of next 
jump time agrees with the one in the fine grid. To do this let us first introduce 
the change of measure for Poisson processes (see section 9.3 of |CT04j ). 

Suppose N t ~ Poi(Ai) under some probability density Pi, then under 
probability density P 2 defined by 

HP A 

^ = exp((A 1 -A 2 )t)(^ (15) 

we will have N t ~ Poi(Aa), where the above term is called Radon-Nikodym 
derivative. 
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In the context of multilevel approach under state-dependent intensity 
model with jump-adapted scheme, we want to calculate 

E[P/]-E[Pti], 

which we can rewrite as 



We shall explain the meaning of this formula: instead of defining in the 
jump- adapted grid formed by cumulative intensity approximated by coarse 
timestep 2he, Pe-\ is defined in the jump-adapted grid formed by fine timestep 
he approximation. 

The Radon-Nikodym derivative Re is defined by 

^ HP . \c 

Re{T) = ^ = exp(Aj - A|) J] ^, 

dF f k A k 

in which T is maturity, denoting \{ = \(S? , U), = X(Sf, tj), A{ = Y^iLi 
and A1 = Yli=i K^i are approximated cumulative intensities upto the ma- 
turity in the fine and coarse grid, k is the index that t[ +1 is the jump time, 
respectively 

This construction is valid since EfP^Lj] = E,[P/_ 1 Re], which can be seen as 
a piecewise constant extension of ffToT) . This can be justified by the theorem 
2.31 in Chapter 2 of |Kar91j . Detail will be done in the future work. 



5.1.2 Variance convergence order 



In the following we shall give some intuitive analysis of the variance conver- 
gence order, which is not a rigorous proof. Hopefully we can use extreme 
path theory to prove it thoroughly in the future work. 

The variance of estimator P# — P^R^ is 



V[P e - Pe-iRe 



< 



Y[Pe-Pe-i + Pe-i(l-Ri)} 



Y[P e -Pe-i] 2 + (WiPe-til-Re)} 



The first part has the order of 0(h). To see the order of Re, let us examine 
the asymptomatic order of squared variance of two components of Re: the 
exponential part and product part alternatively: 
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where X{ = \(S{ , U), XI = X(Sf, t t ). 

We will concentrate on the effect caused in midpoint interval. 



Let h be the timestep of uniform fine grid. Given a midpoint interval 
[tj,tj + 2], since by definition A^ +1 = A^, the exponential part for [i,-, t,- +2 ] is 



j'+i 



E^-^ = (^-^ c j )(h j + h j+1 ) + (Xj +1 -Xj)h j+1 



dX 



0(h)h + —(S(t j )f,t j )(Sf +1 -Sf)h 
O(0 2 ). 



The first term is implied by 0(h) strong convergence and second one 
comes from Sj +1 - Sj ~ SjAW hj+1 ~ 0(h 1 ' 2 ). 

Summing all intervals in the grid, since there are N = T/2h midpoint in- 
tervals, and the rest contributes higher order, asymptotically we get Y^iL\(^{ ~ 
Xf)hi ~ 0(h). Thus the exponential part satisfies 



exp 



Tit 



1 + O(h). 



For the product part, when t,- +2 is a jump time, by similar argument it 
holds that : 



A j+1 Aj- 



l + 0(h 1/2 ). 



So the product part satisfies 

A 
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Asymptotically we get 



V 



Re 



0(h). 



Therefore we have evidence to support that the variance convergence 
order for the multilevel estimator should be 



V[Pe-P e -iRe]=0(h). 

Although cumulative intensity method can deal with the case where jump 
rate is not bounded, the variance convergence order 0(h) leads to a compu- 
tational complexity of 0(e~ 2 (loge) 2 ). 



5.2 Thinning method 

The idea of the thinning method is to construct a Poisson process with a 
constant rate A sup which is an upper bound of the state-dependent rate. 
This gives a set of candidate jump times, and these are then selected as true 
jump times with probability \(S t ,t)/ A sup . 



5.2.1 Algorithm 



Suppose that we have simulated the jump time grid J = {ri, t?, . . . , r m } 
generated by a Poisson process with constant rate A sup , which includes times 
at which jumps occur in [0, T]. On the other hand, consider a fixed time 
grid constituted of N timesteps, t\ — % x ^, % — 1, . . . , A, which is used in 
discretisation schemes for diffusive SDEs. Now the superposition of them will 
be a jump-adapted thinning grid T = {0 = t < t\ < t 2 < . . . < t M = T}. 

For a process which we can simulate the exact increments we have the 
following thinning procedure: 



1. Generate the waiting time for next jump time Tj+i — Tj from a Poisson 
process with constant rate A sup ; 

2. Simulate the evolution of the process up to time T i+ i, 

3. Draw a uniform random number U ~ [0, 1], 
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(a) If p = — - > U, accept r i+1 as a real jump time and 

^sup 

simulate the jump; otherwise go to 2. 

For the general processes we use certain discretisation scheme, e.g. we 
have the following jump- adapted thinning (Milstein) scheme: 

1. Generate the jump-adapted time grid for a Poisson process with con- 
stant rate A sup ; 

2. Simulate each timestep using the (Milstein) discretisation; 



3. When the endpoint t n+ \ is a candidate jump time, generate a uniform 

\(S(t—) r) 

random number U ~ [0, 1], and if U < p T = — , then accept 

^sup 

t n+ i as a real jump time and simulate the jump. 



If we look from the perspective of random measure notation, the thinning 
method can also be formulated in the following ways. First let us recall some 
definitions of point processes and random measures. 

5.2.2 Point processes and random measures 

There are two kinds of basic stochastic processes. The path evolution of the 
first kind is driven by continuous increment, while in the second kind the 
path will only change in the certain jump times. In order to describe the 
stochastic processes of discrete paths, we need to introduce point process. 

Given a filtered probability space (Q,J- t ,F), if we have a increasing se- 
quence of increasing stopping times 

= T < 7\ < T 2 < . . . 

and 

lim T n = oo, 

n->oo 

then the point (counting) process N t associated with stopping times 
(jump times) is defined as 



N t = yi i {T n <t}- 



n>l 
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It counts the number of jumps up to time t. To depict the jump ampli- 
tude (mark) at each jump time, we can define marked point processes and 
associated random measure. 

Suppose mark Y n are random variables taking values on a mark space 
E C R r \ {0}, and Y n are T? n measurable, then (Y n ,T n ) is called a marked 
point process on fix [0, oo]. 

For any A C B(E) and any wefi. 



fi(u; [0,t),A) := ^1{T„( W )<*}1 



{Y„MeA} 



n>l 



defines the random measure associated with marked point process (Y n , T n ). 
It counts the number of jumps within [0, t) whose amplitude belonging to A. 
The class of marked point process is quite general, actually it is a bijection 
to the class of cadlag process. 

For each u e ft, \i (to; -, •) is an Radon measure on (E x [0, oo], B (E x [0, oo])) 
Thus for each P-measuable function /, the integral 

/ f(uj;z,s)/i(u};dz,ds) := / (Y n (u) , T n (u)) 

T„(u))<t 

is a well-defined random variable. 

The compensator ip(dz)ds of the random measure is defined so that for 
all bounded P-measurable function /, we have 

/ f(z, s)fi(dz, ds) - / / f(z,s)ip(dz)ds 

Jz£E JO J z£E 

is a martingale with respect to Tf 

Following the assumption of Platen, we assume that j z£E <f(dz) < oo for 
all s > 0. 

Under those notations, (T2J) can be rewritten as: 



dS(t) = a(S(t-),t)dt+b(S(t-),t)dW(t)+c(S(t-),t) / (z-l)p x (dz, dt), < t < T. 

Jz&E 
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where p\ (u; -, ■) is a Poisson random measure, which means the compensator 
tp(dz) = Xg(z)dz is time independent and g(z) is p.d.f of the mark . In this 
case the waiting time T n+ \ — T n is exponential distributed with a parameter 
A. 

Those general definitions are to allow more flexibility of the dynamics 
of SDEs, e.g. it can admit state-dependent intensity. The dynamics of the 
state-dependent jump-diffusion SDEs we will deal with can be written as 

dS(t) = a{S{t-),t)dt+b{S{t-),t)dW{t)+ [ c{S{t-),t, z)f*{dz, dt), < t < T. 

(16) 

The compensator of /i is defined to be ip(S(t— ), dz)X(S(t— ), t)g(z)dz. We 
also adopt the assumption in [?] that it is bounded by a constant A sup and 
is absolutely continuous. 

The jump-adapted thinning Milstein scheme for (TIB"]) can be formulated 

as 



S 



n+1 



S n + a n h n + b n AW n + \ b' n b n (AW* - h 



S 



n+1 



s 



n+1 



G£ { 



A(S" 



.'n + l) 



>Ui} 



c(S n+1 ,t n+ i, z)fi(dz, t n+ i). 



5.2.3 Multilevel treatment 

In the multilevel implementation, if we use the above algorithm with different 
acceptance probabilities for fine and coarse level, there may be some samples 
in which a jump candidate is accepted for the fine path, but not for the 
coarse path, or vice versa. Because of first order strong convergence, the 
difference in acceptance probabilities will be 0(h), and hence there is an 
0(h) probability of coarse and fine paths differing in accepting candidate 
jumps. Such differences will give an 0(1) difference in the payoff value, and 
hence the multilevel variance will be 0(h). A more detailed analysis of this 
is given in [XGllj . 

To improve the variance convergence rate, we use a change of measure so 
that the acceptance probability is the same for both fine and coarse paths. 
This is achieved by taking the expectation with respect to a new measure Q: 

E[p e - = E Q [p t uri- h-i n K\ 

T T 
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Figure 8: European call option with path-dependent Poisson rate using thin- 
ning without a change of measure 



where r are the jump times. The acceptance probability for a candidate jump 
under the measure Q is defined to be | for both coarse and fine paths, instead 
of p T = X(S(t—),t) / A sup . The corresponding Radon-Nikodym derivatives 
are 

\H, '*U<\; (2p c T , itU<l; 

Rl=l 1 R r={ 1 

[2(1-^), if C/ > - , (2(1- P C T ), if I7> - , 

Since EJ T — R C T = 0(h) and Pi — Pe-i = 0(h), this results in the multilevel 
correction variance Vq[P^ Yl T ~ Pt-i Yl T ^t) being 0(h 2 ). 



The weak convergence of the jump-adapted discretisation with the thin- 
ning procedure is proved in |GM04j , for a class of payoffs on which they im- 
pose to be 4th different iable and that its up to 4th order partial derivatives 
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Figure 9: European call option with path-dependent Poisson rate using thin- 
ning with a change of measure 



are uniformly bounded. By Stone-Weierstrass theorem we can construct a 
sequence of smoothing polynomials which uniformly converges to continuous 
payoff. The limits of this approach is that it is invalid for discontinuous 
payoffs. To prove the convergent variance of multilevel estimator, assuming 
the Lipschitz condition on A, we decompose the estimator into the constant 
rate part and the Randon-Nikdym derivative part, disentangling the effect of 
path-dependence of intensity from the estimator. Under such assumption, we 
can obtain the weak convergence of the estimators for various payoffs by the 
same decomposition, circumventing the difficulties caused by discontinuous 
payoffs. The advantage of this argument is that it can reduce the analysis to 
the constant rate case so that the proof is simplified. 

If the analytic formulation is expressed using the same thinning and 
change of measure, the weak error can be decomposed into two terms as 
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Figure 10: European call option with path-dependent Poisson rate using 
cumulative intensity method 



follows: 



E 



Q 



p,\[r{-p\[r t 



+ E 



Q 



Using Holder's inequality, the bound max(i? r , R{_) < 2 and standard results 
for a Poisson process, the first term can be bounded using weak convergence 
results for the constant rate process, and the second term can be bounded 
using the corresponding strong convergence results |XG11] . This guarantees 
that the multilevel procedure does converge to the correct value. 



5.2.4 Numerical results 

We show numerical results for a European call option using the underlying 
dynamics under risk- neutral measure: 
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dS{t) 
S(t-) 



rdt + adW(t) + + zfi(dz,dt)- zf(z)dz\dt, < t < T, 

J zGE JzeE 



where the random measure \i has compensator \dt with A = 1+ ( S (t-)/s ) 2 • 
The mark has a log normal distribution the density function of which is 
denoted by f(z). 

We use A sup = 1 to generate thinning process. All other parameters as 
used previously for the constant rate cases. 

Comparing Figures M and M we see that the variance convergence rate is 
significantly improved by the change of measure, but there is little change in 
the computational cost. This is due to the main computational effort being 
on the coarsest level, which suggests using quasi-Monte Carlo on that level 
|GW09j . 

The bottom left plot in Figure [H] shows a slightly erratic behaviour. This 
is because the 0(hp) variance is due to a small fraction of the paths having an 
0(1) value for Pf—P^i. In the numerical procedure, the variance is estimated 
using an initial sample of 100 paths. When the variance is dominated by a 
few outliers, this sample size is not sufficient to provide an accurate estimate, 
leading to this variability. 

For comparison, we also show the numerical result given by cumulative 
intensity method in Figure CEDJ The bottom left plot indicates the 0(he) 
variance and the rest plots can be understood consequently. 



6 Conclusions 

In this work we extend the Multilevel approach to scalar jump-diffusion SDEs 
using jump-adapted schemes. The second order variance convergence is main- 
tained in the constant rate case, by constructing estimators using a previous 
Brownian interpolation technique. In the state-dependent rate case, we use 
thinning with a change of measure to avoid the asynchronous jumps in the 
fine and coarse levels. We have also investigated an alternative approach 
which can handle cases in which there is no upper bound on the jump rate. 

The first natural future work is to do rigorous numerical analysis on the 
convergence of variance of correction terms, and weak convergence of the 
ML estimators, which is work in progress |XGllj . The second direction is 
to investigate other cases of model based on specific infinite activity Levy 
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processes, e.g. variance gamma. We also plan to investigate whether the 
multilevel quasi-Monte Carlo method will further reduce the cost. 
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